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We consider optimization of phase response curves for stochastic synchronization of non- 
interacting limit-cycle oscillators by common Poisson impulsive signals. The optimal functional 
shape for sufficiently weak signals is sinusoidal, but can differ for stronger signals. By solving the 
■ Euler-Lagrange equation associated with the minimization of the Lyapunov exponent characterizing 

synchronization efficiency, the optimal phase response curve is obtained. We show that the optimal 
fH , shape mutates from a sinusoid to a sawtooth as the constraint on its squared amplitude is varied. 

I—* 



I. INTRODUCTION 

Synchronization of non-interacting rhythmic elements by common random driving signals [H4l5j , termed stochastic 
or noise-induced synchrony, may explain synchronous behavior of various systems ranging from lasers [|| and electronic 
circuits 0, flll j to spiking neurons [l|, Q and ecological populations [f| , where direct mutual interaction among the 
elements does not exist or is not appropriate to assume. Recent studies have revealed that such synchronization 
generically occurs in a wide class of rhythmic systems including limit-cycle, chaotic, or stochastic oscillators, and also 
for various types of stochastic signals such as Gaussian, Poisson, and chaotic noise 
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Efficiency of the stochastic synchronization is usually quantified by the Lyapunov exponent averaged over noise, 
which measures the mean exponential growth (or decay) rate of small differences between the oscillator states subjected 
to the common noise. For a limit-cycle oscillator und ergoing regular periodic oscillations, the Lyapunov exponent can 
be calculated from the phase response curve (PRC) [18l - l22| . which is a fundamental quantity that characterizes the 



' oscillator dynamics and which has been measured experimentally in many rhythmic elements [23W27] . It can be shown 
fS| . that for limit-cycle oscillators driven by weak Gaussian or Poisson noise, the Lyapunov exponent is always negative 
irrespective of the precise shape of the PRC, ensuring that synchronization always takes place 

What PRC shape yields the best synchronization? For weak Gaussian driving noise, Abouzeid and Ermentrout (28j 
\Q • obtained the optimal functional shape of the PRC by minimizing the Lyapunov exponent with constraints on its 
amplitude and smoothness, which was nearly sinusoidal with a pair of positive and negative lobes (called Type-II, a 
normal form of the PRC near the Hopf bifurcation (2l|). However, the optimal functional shapes may differ for other 
driving signals. 

Here we consider Poisson random impulsive signals with low frequency, which can also induce synchronization 
of limit cycles [Iol - [l3j . When the intensity of the impulse is weak, the linear Gaussian approximation holds and 
the optimal PRC can be shown to be sinusoidal, but for stronger impulses, the optimal solution may take different 
shapes. By using the shooting method [2!| to numerically solve the Euler-Lagrange equation [3(| associated with the 
minimization of the Lyapunov exponent, we show that the optimal PRC gradually deviates from the sinusoid and 
approaches a sawtooth as the constraint on its squared amplitude is varied. Correspondingly, the Lyapunov exponent 
becomes more negative and tends to diverge. Our result implies the importance of nonlincarity in the phase response 
and may provide insights into real- world oscillators such as spiking neurons. 

This article is organized as follows: In Sec. II, some basic facts on synchronization of limit cycles by common 
Poisson noise are presented. In Sec. Ill, we solve the optimization problem and show the gradual transition of 
the optimal solution between sinusoidal and sawtoothed shapes. Section IV summarizes the article with discussions 
on possible relevance of the results to phase response curves of spiking neurons. The appendix gives details on 
wavenumber, symmetry, and phase-plane behavior of the optimal solutions. We also show the optimal PRCs for 
stochastic desynchronzation. 
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II. SYNCHRONIZATION BY COMMON POISSON NOISE 



A. Poisson-driven oscillators 



A pair of non-interactin g; id entical limit-cycle oscillators driven by common Poisson noise can be described by the 
following phase model [TolTilj : 

N(t) 

Oiit) = u+ Y^G{e 1 ,c n )S(t-t n ), 

n=l 
N(t) 

e 2 (t) = u; + J2 G ( e 2,c n )5(t-t n ), (i) 

n=l 

under the assumption that the inter-impulse intervals are sufficiently large such that the oscillator orbit perturbed 
by an impulse relaxes back to the original limit cycle before receiving the next impulse. Here, #12 € [0, 1) are phase 
variables of the oscillators, ui is their natural frequency, N(t) is a Poisson process of rate A, {ii, £2, ■ ■ ■ } are arrival times 
of the Poisson impulses, {c\, c 2 , ■ ■ ■ } are intensities of the impulses (including negative values representing opposite 
directions) independently drawn from an identical probability density function (PDF) P(c), and G(6,c) is the PRC 
of the oscillators. 

The PRC G(9, c) quantifies the asymptotic phase difference of the orbit that is perturbed at phase 9 by an impulse 
of intensity c from the unperturbed orbit [12]. We assume that the PRC G(9,c) is a sufficiently smooth function 
with continuous derivatives G'(9, c) — dG(9, c)/d0, G"(9, c) = d 2 G(9, c)/d9 2 , ■ • • , all of which are periodic in 9, i.e., 
G(9 + Lc) = G(9), G'(9 + 1, c) = G'(8,c), ■ ■ ■ . Equation (QJ is stochastic and should be interpreted in the Ito 
sense (HI- Namely, on arrival of an impulse at phase 9, the phase discontinuously jumps from 9 to 9 + G(9, c) (Tlj |. 



B. Lyapunov exponent 

In Refs. [ljj, EJ, the phase equation (P) is derived from general limit-cycle models by the phase reduction method p~8l — 
l20l |. The Lyapunov exponent A, which quantifies the exponential growth rate of small phase differences between the 
oscillators A9(t) = 9i(t) — Ozit), is given in terms of the PRC as 

A = X J o d9P{9) J dcP(c) In |1 + G'(9,c) \ , (2) 

where P{9) is a stationary PDF of the phase 9 given by a stationary solution of the Frobenius-Pcrron equation 
corresponding to Eq. §T§ jlol - fl2l ] . The phase difference |A0(t)| grows as |A0(i)| ~ |A0(O)| exp(At) when it is small, so 
that the two oscillators tend to synchronize if the Lyapunov exponent A is negative. 

We assume that the impulses are sparse, i.e., the Poisson rate A is small. It can then be shown that the stationary 
PDF of the phase 9 can be approximated as P{&) = 1 + Q(X/ui), so that we may put P(9) = 1 when A is small enough. 
Thus, the Lyapunov exponent is approximately given by [101 lll| 

A = A^ dff J dcP(c) In 1 1 + G'(9,c)\ . (3) 

Moreover, for a sufficiently smooth PRC satisfying G'(9, c) > —1, Eq. © can be bounded from above as 

A <A / dcP(c) f d9G , (9,c) (4) 



=A J dcP(c)[G(l, c) - G(0, c)] = (5) 

by using the inequality ln(l +x) < x and the periodicity of the PRC, so that A is always negative (equality holds only 
for non-physical constant PRCs). Thus, the two oscillators subjected to weak common Poisson noise always tend to 
synchronize. Hereafter, we try to find the optimal PRC that gives the most negative Lyapunov exponent. 

Note that when A is not sufficiently small, we may consider perturbation expansion of the stationary PDF from 
the uniform distribution like P{9) = 1+ (\/ui)P-i (9) + (X/lo) 2 P 2 (9) + ■ ■ ■ to calculate higher-order corrections for the 
Lyapunov exponent, as performed in (Tol. [ill. |28| . For simplicity, we focus only on the case with sufficiently small A 
in the present study. 
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C. Linear Gaussian approximation 

When the derivative of the PRC G'(9, c) is sufficiently small, we may expand Eq. (|3|) as 

A = \j\e J dcP(c) (g'(6,c)- G '( 6 ^ 2 + ...^ 

-~\J dcP ( c )J Q d9G'(9,c) 2 , (6) 

where we used the periodicity of the PRC. Also, if the impulse intensity c is sufficiently weak, the PRC G(9,c) can 
be linearly approximated by using the phase sensitivity function Z{9), which gives the linear response coefficient of 
the phase to infinitesimal perturbations [l8l - |20j . as 



so that A can be approximated as 



G(0,c) = cZ(6), (7) 



A^-^p- [ ddZ'(9f (<0), (8) 



where (c 2 ) = / P(c)c 2 dc. This expression coincides with the Lyapunov exponent of the phase oscillators driven by 
weak common Gaussian- white noise 6-8]. The same result can also be derived more rigorously as the diffusion limit of 
the Poisson noise in which the impulses tend to be weak (c — > 0) and frequent (A — > oo) while keeping A(c 2 ) constant 
and small (TTJ . In this diffusion limit, P{9) can be approximated as P{9) = 1 + O (A(c 2 )/w) [28|. Therefore, fixing 
the average impulse intensity small such that A(c 2 ) <Cwis satisfied, P{9) ~ I holds even for high-frequency impulses 
with large A. 

III. OPTIMAL PHASE RESPONSE CURVES 

A. Euler-Lagrange equation 

The Lyapunov exponent A is a functional of the PRC G{9, c) or the phase sensitivity function Z{9) as given in 
Eq. ([3|) or Eq. ([8]). We try to obtain the optimal shape of G(#,c) or Z(9) for synchronization by minimizing A with 
appropriate constraints. Let us omit the dependence of the PRC G(9, c) on c for the moment. We try to find the 
minimum of the action [3^] 

S[G] =A[G] + fiJ[G] + vK[G] 

= { L(G(6),G'(9),G"(9))d9, (9) 



where A[G] is the Lyapunov exponent, J[G] and K[G] are two independent constraints on the PRC and its derivatives 
(we consider up to the 2nd order), /i and v are Lagrange multipliers, and L(G(9),G'(9),G"(9)) is a Lagrangian. The 
corresponding Euler-Lagrange equation is given by 

d2 dL d dL dL 0, (10) 



d9 2 dG" d9dG' dG 

where the periodicity of the PRC and the derivative, G(9 + 1) = G(9) and G'(9 + 1) = G'(9), are used to eliminate 
the surface terms. When we consider optimization of Eq. ([5]), the PRC G(9) in the above equations is replaced by 
the phase sensitivity function Z(9). 

B. Linear Gaussian approximation 

We here briefly explain the optimal Z(9) under linear approximation. See Abouzeid and Ermentrout [28| for a 
detailed analysis with various constraints. From Eq. ([8]), the Lyapunov exponent of the oscillator in this case is given 
by 



a [z] = ~ j^dez 1 



{Of, (11) 
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where D = A(c 2 ) corresponds to the intensity or variance of the driving noise. We calculate the optimal shape of Z(0) 
by minimizing A [.Z] under the following constraints: 

J [Z] = f Z(9) 2 d9-Bo = 0, (12) 
Jo 

K Q [Z]= [ Z"{6fd0 - Co = 0. (13) 
Jo 

The first constraint Eq. (fl~2|) fixes the squared amplitude of Z(9) to be Bo, which excludes the possibility of non-physical 
divergent Z{9) yielding arbitrarily negative Lyapunov exponents. The second constraint Eq. (|13[) with parameter Co 
restricts the overall smoothness of Z(9). In most realistic finite-dimensional limit-cycle oscillators, the first Fourier 
mode dominates the phase sensitivity function Z{9), reflecting the circular geometry of the limit cycle orbit in the 
phase space. We thus introduce the constraint Eq. (fT3|) to avoid rapid oscillations and choose physically natural 
PRCs, similarly to Ref. [H. 

Introducing Lagrange multipliers /io and vq, the action to be minimized is given by 

S [Z] =Ao[Z] + voMZ] + voKo\Z\ 



de 

L {Z(6),Z'{6),Z"(6))d6. (14) 

The Euler-Lagrange equation determining the optimal Z{&) is given by 

2v Q Z {i) (9) + DZ"{6) + 2^Z{6) = 0, (15) 

where Z^ denotes the 4th derivative of Z. When fi a > and Uq > 0, we obtain a general solution that satisfies the 
periodic boundary condition Z{9) = Z(0 + 1) as 



Z(6) = a^\ X j^ D l- 16 ^° e + P), (16) 
where a and /3 are constants. Due to the periodicity Z{6) = Z(9 + 1), the coefficient of 9 should be quantized as 



J D ± y/^ 2 - 16^Q 

V 4^ = 2 ™' (17) 

where n is an integer number. The constant a is determined from the first constraint Eq. (|12[) as 

1 Z(9) 2 d9 = ^ = B , (18) 
o 1 

namely, a = \fTB~o~. The constant j3 is determined from the boundary conditions for Z{9). Without losing generality, 
we can assume that Z(0) = and Z'(0) > 0, which yields f3 = 0. The second constraint Eq. (fl5|) gives 



Co- (19) 



Equations (|17j) and (|19j) give the relation between Lagrange multipliers (/xq, vq) and the parameters {Bo, Co). In the 
following, we will control the Lagrange multipliers to find optimal solutions with given squared amplitude and overall 
smoothness. 

The optimal phase sensitivity function is thus given by 

Z{9) = v / 2B^sin(27r?i6'), (20) 

which is always sinusoidal regardless of the constraint parameters. The corresponding Lyapunov exponent is obtained 
from Eq. (jHJl as 

DBn n 

Ao = ^n 2 , (21) 
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which decreases with the wavenumber n without bounds. Namely, rapidly oscillating Z(9) can yield very small Ao if 
the constraint on the smoothness of Z(9) does not exist. The second constraint Eq. (TIB")) restricts the range of the 
wavenumber n. In particular, when vq is sufficiently large, only small n is allowed (See appendix). To obtain realistic 
PRCs, we thus set the parameter vq > large enough and focus on Z(9) with n = 1 as well as the corresponding 
G(9), namely, we look for the optimal PRC having only a single pair of positive and negative lobes (Type-II) that 
oscillates only once in 9 <G [0, 1) and crosses the #-axis exactly twice, which is typical of realistic limit-cycle oscillators. 



C. Poisson impulses 



What is the optimal shape of the PRC when the oscillators are driven by common Poisson noise? As we saw, when 
the applied impulse is sufficiently weak and the amplitude of the PRC is small enough, linear Gaussian approximation 
holds and the optimal PRC is sinusoidal. But linear approximation may not be valid when the impulse intensity 
is increased (Io| . On the other hand, if no constraint is imposed on the PRC, an obvious optimal solution is a 
sawtooth, consisting of a straight line of slope —1 and a sharp jump to satisfy the periodic boundary conditions. The 
corresponding Lyapunov exponent diverges to — oo, because a single impulse can already synchronize the oscillators 
by instantaneously rescting their phases to the same value. However, if the impulse is not sufficiently strong to kick 
the oscillator, such a simple solution is impossible. How does the optimal PRC behave in between the two limiting 
situations? 

In the following, we focus on two simple cases in which the oscillators are driven by (i) excitatory impulses with 
a constant intensity (all impulses take the same intensity c), and (ii) both excitatory and inhibitory impulses (the 
impulses take either c = a or c = — a with equal probability). We examine how the optimal PRC deviates from the 
sinusoid and eventually approaches the trivial sawtooth shape as the constraint on the squared amplitude of the PRC 
is increased. 



1. Excitatory impulses 

We assume that the impulse intensity c always takes the same value and simply denote the PRC corresponding to 
this value as G{9). The Lyapunov exponent is given by 

A 1 [G}=\ [ ln\l + G'{6)\d9. (22) 
Jo 

We minimize Ai[G] under the constraints on squared amplitude and overall smoothness of G, 

J[G] = f G{9) 2 d9 - B = 0, (23) 
Jo 

K[G] = ( G"{9fd9 - C = 0, (24) 
Jo 

and examine the dependence of the optimal PRC on the parameter B that determines the squared amplitude while 
fixing C small enough (actually taking the value of v appropriately large) such that the PRC keeps a given level of 
smoothness. 

Introducing Lagrange multipliers fi and v, the action to be minimized is given as 

5i[G] =Ai[G] + fiJ[G] + vK[G] 

= f {A In |1 + G'(9)\ +iu(G{9) 2 - B) + v(G"{9) 2 - C)} d9 
Jo 

= [ Li{G,G',G")d6. (25) 
Jo 

The optimal PRC G(9) is determined by the Euler-Lagrange equation 

d 2 din. d dL x 8L X 



d9 2 dG" d9 8G' dG 



= 0, (26) 
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which gives 

" g(4) + ^(TT^F + mG = ' (27) 

where G*- 4 - 1 denotes the fourth derivative of G. 

If the squared amplitude of the PRC B is sufficiently small, linear approximation for the PRC should hold, i.e., 
G{9) = eZ(9) where e (oc VB) is a small constant. The constraints Eqs. (|2"3")l and become equivalent to Eqs. (TT2")) 
and (|13|) under the linear approximation by rcscaling the multipliers as /i = ug/e 2 and f = vo/e 2 . Substituting these 
into Eq. (j2"T)) . we obtain 

l/oZ(4) + ^-(TlS^ + ^ z = ' (28) 

and taking the e — > limit with D = Ae 2 fixed, we obtain the Euler-Lagrange equation (|T5|) for weak Gaussian noise 
and thus yields sinusoidal Z(9) and G(9) as the optimal solution. On the other hand, if we ignore the constraint 
Eq. (f2"3")) . G(9) = — 9 + const, is a trivial solution to Eq. ([27]) . which gives a sawtooth. Thus, when the squared 
amplitude of G(9) is controlled, mutation of the optimal PRC between the two limiting shapes is expected. 

To confirm this, we numerically calculate a family of solutions to Eq. ([27]) using the shooting method (29[. Namely, 
we numerically integrate Eq. (|27|) by the Runge-Kutta method with adaptive time grids from 9 = to 9 = 1 and 
find appropriate initial conditions G(0), G"(0), G"(0) and G'"(0) satisfying the periodic boundary conditions at 
9 = and 9 = 1. We vary the Lagrange multiplier /i > 0, obtain the corresponding optimal PRC, and check if its 
squared amplitude was equal to the constraint B. Solutions to Eq. (f27|) exist also for a < 0, but they maximize the 
Lyapunov exponent rather than minimize it, and thus are optimal not for synchronization but for desynchronization 
(see Appendix). It can be shown that large values of v lead to small wavenumber (long wavelength) solutions (see 
Appendix). We fix the multiplier v at v = 10~ 5 , which is large enough, to choose non-trivial solutions that cross 
the 9-axis exactly twice in [0, 1) corresponding to the n = 1 case in Eq. (|2"D)) . No periodic solutions are found when 
v < 0. Properties of the optimal solution can be well understood by approximate phase-plane analysis as explained 
in Appendix. 

Figure QJa) shows the results, where the optimal solutions are wrapped within the range [—0.5, 0.5) by taking 
modulo 1. All solutions lay within the plotted region, and no other solutions outside of this region are found. The 
solutions are symmetric with respect to = 0.5 reflecting the symmetry of the Euler-Lagrange equation (f2~T)) (see 
Appendix). As expected, we see that the optimal PRC is almost sinusoidal when the parameter B is small. As B is 
increased, the optimal PRC gradually deviates from the sinusoid and approaches a symmetric sawtooth limit (which 
gives B = 1/12). Correspondingly, the Lyapunov exponent Ai plotted in Fig.[TJb) becomes more negative and tends 
to diverge, and its inverse t% = — 1/Ai, which gives characteristic time for the stochastic synchronization, gradually 
decreases to zero as shown in Fig. []Tc) . 

The optimality of the obtained PRC can be clearly demonstrated by numerical simulations. Figure [2] shows real- 
izations of the stochastic synchronization processes with the optimal and suboptimal (sinusoidal) PRCs. We see that 
the stochastic synchronization occurs much faster when the optimal PRC is used. 

2. Excitatory and inhibitory impulses 

We next consider the case that the intensity of the impulses takes two values ±a with equal probability, namely, 
P(c) = [5(c — a) + S(c + a)] /2. The Lyapunov exponent is 

A=A f ^-{5{c-a) + 5{c + a)}ln\l + G / (9 7 c)\d9dc 
Jo 2 

=^J \n\(l + G' (9, a)) (1 + G' (9, -a))\d9. (29) 

For simplicity, we seek for symmetric PRCs that satisfy G(9, —a) = —G(9, a). This condition should be always satisfied 
if a is sufficiently small, because the PRC can be linearly approximated as G(9,c) = cZ(9). The existence of the 
diffusion limit is also ensured with this condition [Tl|. Note that, for stronger impulses, the PRC generally becomes 
asymmetric and does not satisfy the above condition. We here focus only on the symmetric case for simplicity. 
The Lyapunov exponent is then given by 

A2[G] = ^ / \n\l-G\9) 2 \d9 (30) 
2 Jo 




FIG. 1: (Color online) (a) Numerical solutions of the Euler-Lagrange equation (|27[) obtained by the shooting method. The 
dashed line plots the limiting sawtooth solution. The solid curves are non-trivial solutions for various values of the squared 
amplitude ranging from B = 2.98 x 1CP 4 to 7.07 x 1CP 2 . (b) Dependence of the Lyapunov exponent Ai on B. (c) Dependence 
of the characteristic synchronization time r\ = — 1/Ai on B. 



with the abbreviation G(9) = G(9, a). We minimize A2[G] under the constraints (f2"3"f and (jM)) . Introducing Lagrange 
multipliers fi and is, we obtain the action 

S 2 [G] =A 2 [G] + fiJ[G] + vK[G] 

A In 1 1 - G'(9) 2 1 + n (G(9) 2 -B)+v (G"{0) 2 - C)\ d9 



f L 2 (G,G',G")d6, (31) 
Jo 



and the associated Euler-Lagrange equation 



A G"(l + G' 2 ) 



(9) + - (1 l _^ /2)2 y + i-iG = 0. (32) 



If the squared amplitude B of the PRC is sufficiently small, we can rewrite Eq. (|32j) using the linear approximation 
of the PRC with rescaled multipliers, G(9) = G(9, a) = aZ(8), fx = Ho/a 2 and v = v$/a 2 , as 
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FIG. 2: (Color online) Numerical realizations of the stochastic synchronization processes with the optimal and suboptimal 
(sinusoidal) PRCs for the case of excitatory impulses. The squared amplitude of both PRCs is set at the same value B = 0.045. 
(a) Optimal PRC (Lyapunov exponent Ai = —1.221). (b) Sinusoidal suboptimal PRC (Lyapunov exponent Ai = —0.025). 
(c),(d) Numerical realizations of the stochastic synchronization processes with the PRCs shown in (a),(b). 



Taking the diffusion limit, i.e., a — > and A — >■ oo with D = Xa 2 fixed, the Euler-Lagrange equation (|15[) under 
the linear Gaussian approximation is derived. Therefore, we obtain a sinusoidal Z(9) and hence G(9) as the optimal 
solution for small B. On the other hand, if we ignore the constraint, Eq. (f52j) has the obvious solution G(9) = —6 as 
before. In the present case, additionally, G(9) = 9 is also an optimal solution because G(9, —a) = — G(9, a). 

Using the numerical shooting method, we obtain a family of optimal solutions to Eq. (|32| as plotted in Fig. [3ja). 
As in the previous case, the multiplier v is fixed at 10 -5 , which is large enough to yield smooth PRCs. Unlike the 
previous case, no solution with period 1 exists when /i < 0. As the parameter B increases, the optimal PRC gradually 
deviates from the sinusoid. In this case, the PRC approaches a double sawtooth, in contrast to the single sawtooth 
that we obtained previously, reflecting the symmetry assumption. The Lyapunov exponent A2 becomes more negative 
and tends to diverge, and the characteristic synchronization time T2 decreases to zero as shown in Figs.[3Jb) and (c). 
The optimality can be demonstrated by numerical simulation as shown in Fig. [¥] 



IV. DISCUSSION 



We considered the optimization problem of the PRC for synchronization of limit-cycles oscillators by common 
Poisson noise and observed a crossover of the optimal PRC from a sinusoid to a sawtooth by increasing its squared 
amplitude. Now we take some time to stress the importance of considering nonlinear PRCs. The phase sensitivity 
function Z(9) quantifies the linear response property of the oscillator phase to infinitesimal perturbations, which is 
determined by the local phase-space structure of the oscillator near the limit-cycle orbit jl8l - i2C| . In contrast, the PRC 
can reflect nonlinear dynamics of the oscillator away from the limit-cycle orbit by finite distances, providing more 
detailed information. Also, in many experiments, applied perturbations to the oscillator arc not always sufficiently 
small and nonlinear effects can become important. In the present study, we considered only two simple types of 
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FIG. 3: (Color online) Numerical solutions of the Euler-Langrange equation (|32|l obtained by the shooting method, (a) The 
dashed line plots the trivial solution. The solid curves are non-trivial solutions for various values of the squared amplitude 
ranging from B = 1.04 x 10 -3 to 1.70 x 10 -2 . (b) Dependence of the Lyapunov exponent A2 on B. (c) Dependence of the the 
characteristic synchronization time T2 = — I/A2 on B. 



driving impulses, i.e., (i) excitatory and (ii) excitatory and inhibitory impulses, and also assumed symmetry of the 
PRCs in the latter case. More general types of driving impulses and asymmetric PRCs can be considered within the 
same framework, though they are beyond the scope of the present study. For example, it would be interesting to seek 
for the optimal family of PRCs for a given distribution of the impulse intensity c by making additional assumptions 
on the c-dependence of the PRC G(9,c). 

Are there examples of the optimal PRC in nature? In neurophysiology, the PRCs of periodically spiking cells 
have been recorded in many experiments p3 - l27l |. For example, Tateno and Robinson [25| calculated the PRCs of 
periodically spiking interneurons from monkey somatosensory cortex and examined their dependence on the intensity 
of applied perturbations. As the intensity increases, the PRC changes its shape from sinusoidal to sawtoothed (Figs. 4 
and 5 in Ref. |25|). This dependence of the PRC on the applied signal intensity resembles the gradual transition 
that we obtained in Fig. [TJ The authors also found that the sawtooth-like PRCs lead to faster synchronization of the 
neurons (2(| . Nesse and Clark (27j calculated the PRC of photoreceptor cells from marine invertebrate Hermissenda 
and revealed noticeable linear dependence of the PRC on 9 (Figure 5 in Ref. [13 )■ The authors suggested that the 
reset effect of such a PRC may be helpful for network information processing. 

Stochastic synchrony can be a mechanism for long-range synchronization of gamma oscillations in the cortex @]. 
Because the thalamus is at the center of the brain and communicates with all cortical regions, it is a good candidate 
to provide common input to areas of the cortex that are far apart and not directly connected. This shared thalamic 
drive represents a straightforward mechanism to mediate synchrony between these areas, and the PRC of neurons in 
cortex receiving thalamic input could be optimized for this purpose. In Ref. [Hj], Lefort et al. report that synaptic 
strength of neurons in some cortical areas have a long-tailed distribution, meaning certain synapses are much stronger 




FIG. 4: (Color online) Stochastic synchronization processes with the optimal and suboptimal PRCs for the case of both 
excitatory and inhibitory impulses. The squared amplitude of both PRCs is set as the same value B — 0.017. (a) Optimal PRC 
(Lyapunov exponent Ai = —0.690). (b) Sinusoidal PRC (Lyapunov exponent Ai = —0.544). (c), (d) Numerical realizations of 
the stochastic synchronization processes with the PRCs shown in (a), (b). 



than others (the amplitude of postsynaptic potentials spans over a few millivolts). Thus, it might actually be more 
appropriate to consider finite-intensity impulses than weak Gaussian noise as the driving signal to the neurons. 

The optimization viewpoint may give interesting insights into the understanding of biological systems, because they 
evolved to perform certain biological functions efficiently. If the stochastic synchronization mechanism is used in some 
biological systems, their PRC may be optimized to best perform synchronization. The sawtoothed PRCs that we 
obtained are not only optimal for the synchronization by common Poisson noise, but they are singular in the sense 
that they lead to instantaneous phase resetting of the oscillators. Thus, it may not be surprising if such a singular 
shape is actually utilized in real biological systems. This parallel between evolutionary optimization and optimization 
for a desired function certainly makes the interpretation of such a singular shape highly suggestive and intriguing. 
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APPENDIX 



In this Appendix, we give detailed discussion on the dependence on the Lagrange multiplier is, symmetry, and 
phase-plane analysis, of the optimal solutions. We also show the optimal PRCs for stochastic desynchronization. 
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A. Dependence of the optimal PRCs on the multipliers 



We find that if the multiplier vq or v is sufficiently large, only small wavenumber (long wavelength) solutions are 
allowed for Z or G. This can be proven for the Euler-Lagrangc equations (|T5|) . (|2"T)) and (|3l?|) . 



1. Linear Gaussian approximation 

We consider a solution of the Euler-Lagrange equation (|15p with wavenumber n and denote the corresponding 
Lagrange multipliers {[i n ,v n ). Substitution into (| 15[) yields 

8Tr 4 v n n 2 + Dn 2 + ^ = 0, (34) 

namely, the Lagrange multipliers scale with the wavenumber n as fi n oc n 2 and z^„ oc 1/n 2 . Thus, larger \i and smaller 
^ lead to PRCs with larger wavenumbers. We set v sufficiently large to obtain the n = 1 solution in the main text. 



2. Poisson impulses 

Rescaling the phase variable as 9 — > n9, the Euler-Lagrange cquaion (|27j) is transformed to 

^ G ^4 {l+ fiZw + ^ G{ne) ^ (35) 

Defining a rescaled PRC G n (9) = G(n9)/n, the above equation can be cast into the same form as Eq. (|27[) . 

^ 4 >(0) + ^ (1+ G ^ ))2 + MnG n (&) = 0, (36) 

where rescaled Lagrange multipliers fi n = n 2 fj, and = v/n 2 are introduced. Thus, if G{9) is a solution of Eq. ([27)1 
with multipliers /i and z^, its rescaled function G n (9) is also a solution of Eq. (|27[) with rescaled multipliers and 
i'ti (n = 1)2, • • •). This implies that larger wavenumber solutions (n > 1) correspond to larger fi and smaller v. As 
shown in Fig. [T] the multiplier /i control the squared amplitude B and determine the shape of the periodic solutions. 
Thus, v determines the wavenumber of the solution, which we take sufficiently large [y = 10~ 5 ) to obtain the PRC 
corresponding to n = 1. Similarly, rescaling Eq. (|32p . we obtain 



(l-G' n (0W 

to find = n 2 \i and = v/n 2 . Thus, if v is sufficiently large, the PRC takes the smallest wavenumber n = 1. 



B. Symmetry of the optimal solution 

From Eqs. (J27|) and © with periodic boundary conditions G(0) = G(l) = 0, G"(0) = G'(l), and G"(0) = G"(l), 
we obtained symmetric solutions with respect to 9 = 0.5 as shown in Figs. QJa) and[3ja). These PRCs inherit the 
symmetry from the Euler-Lagrangc equations (or the actions to be minimized). To see this, let us define a function 
F(9) as 

F{9) = -G{\-9). (37) 

The derivatives are F'{9) = G'(l - 9), F"{9) = -G"(l - 9), and F^(9) = -G< 4 >(1 - 9). Substituting into Eq. (j27| . 
we find that F(9) obeys the same Euler-Lagrange equation as G(9), 

^ (4) W + ^7i-TW + ^- ' ( 38 ) 
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with the same boundary conditions G(0) = F(0), G'(0) = F'(0), and G"(0) = F"(0) (note that all the solutions 
satisfy G"(0) = 0). We thus obtain G(9) = F{9) = — G(l — 9), indicating that G(9) is symmetric with respect to 
9 = 0.5. 

Similarly, when the optimal PRC G(9) obeys the Euler-Lagrange equation (|32p . we can show that F(9) satisfies 
the same equation 

«/fW(g) + ^^ff +MF = (39) 

with the same boundary conditions. Thus, G(9) = F(9) = — G(l — 9) holds and G(9) is also symmetric with respect 
to 9 = 0.5. 



C. Phase-plane analysis 

As we explained, we fix the multiplier v large (but still much smaller than unity, v = 10~ 5 <C 1) to obtain physically 
realistic PRCs. Here, to gain insights into how the shapes of the optimal PRCs are determined, we set v = and 
ignore the 4th-order derivatives in the Euler-Lagrange equations, which does not affect the solutions qualitatively. 
With this approximation, the dependence of the optimal solution on the constraint B or on the Lagrange multiplier 
\x can be clarified by a simple phase-plane analysis. 



1. Excitatory impulses 

We set v = to approximate Eq. (|2"T)) as 



G"{9) = -^G{9){l + G'{9)f (40) 
A 



and rewrite this equation as 



G'{9) =H(6), 

H\9)=-^-G(9){l + H(9)f . (41) 
A 

We examine the orbit of this two-dimensional dynamical system as a function of 9 £ [0, 1) on the G — H plane with 
a periodic boundary condition G(0.5) = G(— 0.5) and 11(0.5) = 0.5). 

Let us assume /i > first. Figure shows an example of the vector field at /i = 10. The horizontal line H = — 1 
is a separatrix corresponding to the sawtooth solution G'(9) = — 1. All orbits starting from H > — 1 are closed, 
implying the existence of a conserved quantity. Applying Noether's theorem [3(| to the Lagrangian in Eq. (|27j) . we 
find that the quantity 

Ci = A { if^7g) - ln|l + G'W|| - fiG(9) 2 (42) 

is actually conserved along the flow generated by Eq. (|4"Tj) , reflecting the translational symmetry of the Lagrangian 
with respect to phase, namely, that the Lagrangian does not depend on 9 explicitly. 

A solution possessing period 1 is chosen from this family of closed orbits by the shooting method. The solid loop 
in Fig. 03a) shows such a periodic solution, and the solid curve in Fig. [SJb) is the corresponding optimal PRC. No 
orbit starting from H < — 1 can form a closed loop, because the vector field points to the upper- right and lower right 
in the third and fourth quadrant, respectively. Thus, in this region, the orbits have to jump from G(9) = —0.5 to 
0.5 as shown by a broken curve in Fig. [S^a). However, the periods of such orbits are always less than 1 and therefore 
solutions with G'(8) < — 1 do not exist. 

It can be seen from Eq. (|41[) that the Lagrangian multiplier /i determines the time scale of the dynamics in the 
vertical H direction. As [i increases, the vertical dynamics becomes faster, so that the orbit is more strongly attracted 
to the separatrix H = — 1 and tends to move along it, as shown in Fig. (5][c). Correspondingly, the optimal PRC 
approaches the sawtooth as shown in Fig. [5jd). Note that the separatrix G' = H = — 1 persists even if v > 0. This 
can be confirmed by taking the limit G' — >• — 1 in Eq. (|2"T)l . which gives G" —> 0. Thus, the sawtooth limit also persists 
in the original system. 

When /i < 0, we obtained the optimal PRCs for desynchronization as summarized in Appendix D. 
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FIG. 5: (Color online) Vector fields and corresponding PRCs for the case of excitatory impulses. Contour color represents 
magnitude of each vector, where darker color corresponds to higher magnitude, (a) Vector field at /i = 10. Solid loop is an 
orbit with period is 1. Broken curve is an orbit starting from H < — 1. (b) Optimal PRCs for /i = 10. (c) Vector field at 
/i = 15. Solid loop plots an orbit whose period is 1, and (d) the corresponding PRC. 



2. Excitatory and inhibitory impulses 

The same analysis can be applied to the case with both excitatory and inhibitory impulses. As shown in Fig. [51 
when ji > 0, horizontal lines H = ±1 are the separatrices. Orbits starting from \H\ < 1 always form closed loops, 
while those starting from \H\ > 1 cannot form a period-1 solution. The conserved quantity in this case is given by 

C-2 = -X L^L 2 + \ In |1 - G'(0) 2 | j - /iG(tf) 2 . (43) 

Increasing the multiplier /x, the optimal solution gradually expands and changes its shape from a circle to a rectangle. 
The corresponding PRC deviates from a sinusoid and approaches a double sawtooth. As before, the separatrices 
persist even if v > 0. No orbit with period 1 was found when /i < 0. 



D. Optimal PRCs for stochastic desynchronization 



The Euler-Lagrange equation gives the solutions that yield the extremum of the action, namely, minimum and 
maximum, of the Lyapunov exponent A under the given constraints. In the case of excitatory impulses, we can 
vary the Lagrange multiplier [i controlling the squared amplitude of the PRC in the negative range, [i < 0, while 
keeping the other Lagrange multiplier v the same as in the main text, v = 10 -5 , to obtain the PRC that maximizes 
the Lyapunov exponent. The corresponding Lyapunov exponent is positive, indicating that the PRC is optimal for 
stochastic desynchronization [ll[. As shown in Fig. [TJa), this optimal PRC has a sharp cusp at 9 = when /j is 
sufficiently negative, and gradually approaches a sawtooth as fi increases. 

Examples of the optimal PRC and the corresponding phase-plane orbit are plotted in Figs. EJb) and[7|c). It is 
interesting to note that the PRC plotted in Fig. [5](b) or (d) is "type 1" while the PRC in Fig. [7ta) is "type 0" in 
Winfree's classification [l9|, HH ; the type 1 PRC is continuous and is observed for moderate perturbation intensity, 
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FIG. 6: (Color online) Vector fields and corresponding PRCs for the case of both excitatory and inhibitory impulses, (a) 
Vector field at fi = 20. Solid loop plots an orbit whose period is 1. (b) Vector field at pi = 100. Solid loop plots an orbit whose 
period is 1. (c) Vector field at /i = —2.5. (d) Optimal PRCs for two values of /i. Solid curve is the result for fi = 20. Broken 
curve is for fi = 100. 



whereas the type PRC is discontinuous and is observed when an oscillator is strongly perturbed [34| . Thus, under 
the present criteria, the optimal PRC for stochastic synchronization is type 1 and that for desynchronization is type 
0. 
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